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ABSTRACT 

The maximum -likelihood sense optimum detector for the detection of randomly 
fluctuating multiple highlight objects in Gaussian reverberation and noise is an Estimator- 
Correlator (EC). If a large time-bandwidth product signal is used, a conventional matched-filter 
detector can spread scattered energy into multiple detector bins, resulting in lower output in each 
bin than if all the energy had been concentrated in just one. 

The wideband, or Wavelet Transform Domain EC (WTD-EC) recombines the detection 
statistics from discrete highlights into a single detection statistic, possibly dramatically improving 
detector performance as characterized by the Receiver Operating Characteristics (ROC). 
Implementation of the estimator branch of the EC requires an a priori estimate of the second 
order statistics of the echo, i.e. the covariance function or scattering function. Scattering 
functions are used as prior statistics because they are physically meaningful, robust statistical 
characterizations of distributed scatterers. A critical issue in the implementation of the WTD-EC 
is the sensitivity of its performance to the errors in the a priori statistical models. 

This thesis analyzes the sensitivity of the WTD-EC to errors in the assumed scattering 
functions. Analytic expressions and verified numerical code are developed for the computation of 
ROC curves for arbitrary target scattering functions, reverberant scattering functions, and 
wideband transmitted signals. The sensitivity of the detector is found by deriving an expression 
for the ROC when the assumed scattering function is not the correct scattering function. 

It is demonstrated that the WTD-EC can not be characterized by a scalar figure of merit, 
and that only a ROC curve is adequate to characterize the performance of the detector. It is 
demonstrated that the sensitivity closely mimics the autoambiguity function of the transmit signal. 
It is shown that, to some extent, one can trade time-bandwidth product for robustness to error. 
Compelling evidence is given that the WTD-EC does in fact retain it’s advantages even when the 
a priori estimate contains large errors. Finally, it is shown that one need only estimate the 
position of target scatterers. It is s not necessary to estimate the reverberation scatterers, or the 
relative magnitude of the target scattering highlights. 
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"Worse than being blind 
is being able to see 
but have no vision." 


- Helen Keller 



Chapter 1 
INTRODUCTION 
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1.1 Motivation 

The process of transmitting a known signal into an environment and listening for its 
reflections is referred to generally as Range-Doppler detection. It is usually known under the 
more specific names SONAR, RADAR, Medical Ultrasound, SODAR, Non-Destructive 
Evaluation, etc. Advances in computer memory, computation speed, and transducer bandwidth 
have made possible (and will continue to make possible) higher Time-Bandwidth product (TB) 
transmit signals and therefore higher resolutions. Traditional Matched Filter (MF) and Estimator- 
Correlator (EC) methods have better performance with higher TB signals, as long as the 
resolution of the signal remains coarser than the size of the target. When the TB becomes high 
enough, resolution becomes finer than the size of the target, the reflected energy is spread across 
several detection bins and can result in lower probability of detection. This is sometimes referred 
to as “splitting loss” or “spreading loss”. These distinctly resolved reflections will be called 
highlights. The purpose of this thesis is to use an EC structure to recombine the multiple discrete 
highlights of the signal, recover the splitting loss that comes with high resolution signals, and 
evaluate the performance of this EC. To do this, an a priori estimation of the target scattering 
function must be made, that is, one must make a guess as to where the highlights are. One of the 
risks in taking on this project is that errors in the scattering function estimate could cancel the 
gains made by recombining the energy, thus making this process useless in practice. Compelling 
evidence will be presented that highlight recombination retains its advantages even when the a 
priori estimate contains large errors. 

1.2 Statement of Problem 

Higher time-bandwidth product signals are desirable because they afford higher resolution 
and better probability of detection under correlated noise. With traditional detectors this is 
limited by splitting loss. A method is needed to recombine the energy split into different detection 
bins by the high resolution of the signal. 

Two detection methods applicable to the Range-Doppler problem are considered. They 
are, the Estimator-Correlator (EC) and the Matched-Filter (MF). The matched filter simply 
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correlates the received signal with replicas of the transmit signal which have been delayed and 
frequency shifted (or scaled). The delay and frequency shift (or scale) of the replica of the 
transmit signal with the largest correlation represents the distance and radial velocity of the target. 
When the change in the reflected signal due to the velocity of the target is modeled as a frequency 
shift, this is called narrowband processing. When the change in the reflected signal due to the 
velocity of the target is modeled as a scaling of the signal, this is called wideband processing. 
Wideband processing is explained in section 2.3, narrowband in section 2.4. 

The estimator-correlator uses the covariances of the received signal when the target is 
present, and when it is absent, to form an estimate of the target signal, which is then correlated 
with the received signal. The EC can also be implemented with either wideband or narrowband 
processing. 

Recently, an implementation of the EC has been derived which represents a spread target 
as a scattering function in the Wavelet Transform Domain (WTD) [1]. Also, Dixon and Sibul 
have done work in parameter estimation for spread targets [2], [3]. In[4] Roan and Sibul derived 
the Receiver Operating Characteristics (ROC) of the WTD-EC, but they only computed the ROC 
for the case where the covariance eigenvalues are all equal, and they did not compute those 
eigenvalues from the scattering functions and autoambiguity function of the transmit signal. It is 
necessary to characterize the performance with the ROC and not with a scalar figure of merit 
because scatterers with different highlight structures but the same total energy have different ROC 
curves, and thus can not be described by a single parameter. Critical issues that remain to be 
addressed include: implementing WTD-EC for arbitrary eigenvalues which have been computed 
from the scattering functions and the autoambiguity function of the transmit signal, numerical 
problems in computing the eigenvalues from the scattering function estimates, and computing the 
Heaviside expansion of the probability distributions, analytic methods of validating the 
computation, computing ROC curves directly from scattering functions and signals of interest. 

One is unlikely to ever estimate the scattering function exactly. It may be true that errors 
in the scattering function estimate cancel the benefits of attempting to recombine the highlights, 
making this process useless in practice. A sensitivity analysis of the WTD-EC needs to be derived 
which should include: general properties of the sensitivity, analytic special cases to validate the 
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sensitivity, and computations of the sensitivity for scattering functions and transmit signals of 
interest. 

1.3 Original Contribution 

The receiver operating characteristics of the wavelet transform domain estimator- 
correlator are computed for arbitrary covariance eigenvalues which have been computed from the 
scattering functions and autoambiguity function of the transmit signal. Numerical problems and 
solutions are identified. Analytic special cases are used to validate the numerical method. 
Performance of the WTD-EC is compared with the traditional (one highlight) estimator 
correlator. It is shown that the WTD-EC does in fact have a gain in performance proportional to 
the energy distributed between the highlights. It is demonstrated that the distance between 
highlights has little effect on the performance. 

The Receiver Operating Characteristics (ROC) are derived for the case where the a priori 
estimate is not correct. This derivation is independent of the wavelet transform domain 
implementation of the EC, and is true regardless of how the eigenvalues of the covariances are 
computed. Anal ytic special cases are derived to validate the numerical implementation. Many 
properties of the sensitivity of the WTD-EC are shown. The sensitivity is related to the width of 
the autoambiguity function. One can trade off, to some extent, time-bandwidth product for 
robustness to error, but there is a minimum TB below which there is no performance advantage. 
The sensitivity does not depend on the input SNR. Compelling evidence is given to show that 
making a scattering function estimate is always better than not making one (assuming one 
highlight). Overestimating the overall magnitude of the target highlights does not significantly 
degrade the ROC curve, and thus, when implementing the detector, one need not make an 
estimate of the reverberant scatterers. Incorrectly estimating the relative magnitudes of highlights 
with respect to each other also does not significantly degrade the probability of detection, unless 
the estimate is so small that it is equivalent to guessing that there is no highlight at all. 


4 


1.4 Assumptions 

It is assumed that the propagation medium is isovelocity, homogeneous, and has no 
frequency dependent attenuation. Scattering is assumed to be wide-sense stationary and 
uncorrelated(WSSUS) [5]. This allows for a convenient matrix formulation connecting the 
wideband scattering function to the covariance matrices. Since no phase information about the 
highlights is necessary, this makes for a robust model for the recombination of the energy split 
into highlights. The wideband scattering model is assumed, which implies that the acceleration of 
the targets is smaller than a limit set by the time bandwidth product, and that the transmit signal is 
finite energy and square integrable. The scattering function model is very comprehensive. The 
target object, channel, and unwanted scatterers are modeled with scattering functions to represent 
highlights, multipath, biologies in underwater, etc. The total scattering function is the convolution 
of all the component scattering functions [6]. In this thesis the target and reverberation are 
modeled with scattering functions, therefore one is not limited to assuming uncorrelated Gaussian 
white noise. 

1.5 Overview 

Chapter 2 reviews the literature on the topics of signal design and wideband sonar 
processing. Prerequisite results and theories in wideband processing and estimation are derived. 

In Chapter 3 various numerical difficulties are overcome, a computational code for calculating the 
ROC curve is implemented, several analytic cases are derived to validate the code, and ROC 
curves are computed and compared to traditional methods. In Chapter 4 the sensitivity of the 
detector to errors in the scattering function estimate is derived, computation code is implemented, 
analytic cases are used to validate it, and the sensitivity of the detector for several interesting 
kinds of errors is computed. Appendix B covers the numerical difficulties of computing the 
Heaviside expansion when poles are close together. 


Chapter 2 
BACKGROUND 
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2.1 Wideband Processing 

One way to model the reflection from a moving target is to say that it imparts a frequency 
shift to the reflected signal. This is called narrowband processing, and is discussed in more detail 
in section 2.4. When the velocity of the target is very high, or when the time-bandwidth product 
of the signal is very large, the amount by which a signal shrinks or stretches in duration (called 
time dilation) becomes significant. A reflection model that includes time dilation is called 
wideband processing, and is discussed in section 2.3. 

Several papers have been written by Weiss showing how wideband processing can be 
implemented as a wavelet transform [7], [8], [9], [10]. Jin et. al. have also shown this [11]. Altes 
has derived a set of signal transforms which leaves the wideband ambiguity function unchanged 
[ 12 ]. 

Jourdain and Henrioux have investigated the wideband auto-ambiguity properties of 
BPSK signals, linear FM, Hyperbolic FM, and Costas Codes [13]. Modulating a BPSK with a 
Gaussian window and repeating small sections of the code produces higher peak to sidelobe levels 
than other signal types. Methods of approximating the wideband ambiguity function have been 
investigated by Ricker [14] and Glisson et. al.[15]. In these schemes scaled versions of the 
transmit signal are pre-fabricated and the received signal is narrow band correlated around each 
replica. 

In “Wavelet Domain Implementation of the Estimator-Correlator and Weighted Wavelet 
Transforms”, Sibul et. al. [1] derive a continuous wavelet domain estimator correlator for signals 
propagated in stochastic scattering channels. They show that the detection statistic must be 
defined using a reproducing kernel Hilbert space. Tague and Sibul [16] have developed a 
canonical Estimator-Correlator array processor. 

In Roan’s dissertation [17] a Kalman filter method is used to track changes in the 
scattering function. ROC curves are not derived, nor is the detection statistic’s sensitivity to error 
in the scattering function estimate computed. 



6 


In Grove’s dissertation [18] the relationship between errors in physical parameter 
measurements and errors in the resulting scattering function estimate are derived. Analysis of the 
change in the detection statistic due to these errors is not conducted in detail. 

Very recently Ricker and Cutezo have published a paper [19] where ROC curves are 
computed for the case where the assumed scattering function is not the correct scattering 
function. They give three examples each with greater than five highlights. This thesis takes the 
approach of looking at examples with a very small number of highlights to clearly show 
relationships such as: errors in position and the autoambiguity function, and errors in highlight 
magnitude and performance. 

In Roan and Sibul’s “Performance Quantification for the Wavelet Transform Domain 
Estimator-Correlator” [4] the ROC is calculated for a high time-bandwidth product signal where 
the noise covariance eigenvalues all have the same value. The calculation of these ROC curves 
did not use an actual transmit signal and target and reverberant scattering functions, but rather 
used eigenvalues that were thought to be representative of a typical system. This thesis expands 
on that work by creating a method for calculating the ROC from the actual transmit signal, target, 
and reverberant scattering functions, and deriving the numerical limitations on that calculation 
imposed by the Heaviside expansion. 

2.2 Signal Design 

Being able to find the sensitivity of the WTD-EC for an arbitrary transmit signal has 
important implications for signal design. Future researchers could use the sensitivity analysis of 
this thesis to find new transmit signals with new qualities. It is important to put such an 
optimization in the context of signals optimized for other qualities. To this end, a review of signal 
design literature is given in Appendix A. 

The results of this thesis show that, to maximize the ability to discern between similar 
scatterers (maximum sensitivity to error), one should use high time-bandwidth product signals 
which have an autoambiguity mainlobe which is as narrow as possible in the direction where you 
expect error. Signals with “thumb-tack” like autoambiguities, such as FSK signals or Costas 
codes, are applicable if one is concerned with errors in both range and velocity. If one is only 
concerned with errors in range, LFM or HFM (velocity tolerant) signals would do well. One 
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needs a signal which has the widest autoambiguity mainlobe for a given time-bandwidth product 
in order to maximize robustness to error. This is not a type of signal which has previously been 
derived, and thus this thesis shows that there may be signals yet to be found which have superior 
robustness to error in the a priori scattering function estimate than any used today. 

The product of the total time-duration of a signal and its total bandwidth, called Time- 
Bandwidth product (TB), is an important characteristic of a signal. For many years, it was a “folk 
theorem” that the number of degrees of freedom of a signal equaled 2TB. This was made 
rigorous by Slepian, Landau, and Pollack [20], [21], [22], [23], [24]. Generally, the higher the 
TB, the sharper the resolution in the time-scale, or time-frequency plane. 

Given a digital signal with N sample points, one could “play it back” through a Digital to 
Analog Converter (DAC) at a low sampling frequency, and get an analog signal with a long time 
duration and low bandwidth. Conversely, one could play it back with a high sampling frequency 
and get a signal with a short time duration and large bandwidth. Because a digital signal with N 
samples has the same TB regardless of the sampling frequency used to play it back through a 
Digital to Analog Converter (DAC), this thesis will usually characterize a signal by its TB, and not 
mention the time duration or bandwidth it would have if transmitted at a particular sampling 
frequency. 

23 Wideband Derivation 

In order to implement an estimator correlator one must have a model for how signals are 
reflected by scatterers. Consider a signal f(t) transmitted from a stationary position and reflected 
from a stationary object. The return signal g(t) will be delayed by a factor x = 2R/c so that g(t) = 
f(t-x). Here R is the distance between the transmitter and reflector, and c is the propagation speed 
of the signal [25], [26], [27]. 

More generally, if R is expressed as a function of time, as in the case of a moving 
scatterer, then the round-trip travel time is: 




( 1 ) 
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Then x(t) can be solved for by differentiating with respect to time (denoted by prime) 


( 2 ) 


In this derivation parenthesis denote arguments of functions and brackets denote groups of 
terms. Assuming the velocity v = R'(t) is nearly constant (i.e., no acceleration) over the length of 
the signal 

2v 


r '(0 = 


c+ v 


(3) 


Integrating and choosing x(x 0 ) = x 0 , yields 


2v 


<0= T 0 + —[t- r 0 ] 


c+v 


(4) 


Substituting into g(t) = f(t - x(t)) yields 


c-v 

g(t) = f(s[t - r]) where s = 


This is the wideband model for signal reflection. It is easy to see that the signal is 
stretched or contracted by the scale factor s. This scale factor is analogous to the scale factor of 
the wavelet transform. It will be shown that wideband correlation is the same as a wavelet 
transform. 

The restriction that the velocity be nearly constant (or nearly no acceleration) is more 
precisely stated as [50]: 


TB 2 < 



(5) 


where T is the time length of the signal, B is the bandwidth, and a is the acceleration of the object. 
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In order to estimate the range and velocity of an object, the received signal g(t) must be 
correlated with hypothesized scaled and time shifted versions of the transmitted signal. The 
maximum of this wideband correlation will correspond to the velocity and displacement of the 
object. The wideband correlation of two functions f(t) and g(t) is given by: 

JVC(s, r) /g = s 2 f g(t)f*(s[t - r])dt, ( 6 > 

where the factor s 1/2 has been added to preserve energy. The wavelet transform is defined as [28] 

w fg (s,T )= 4r (7) 

s s 

It can be seen that the wideband correlation (or wideband matched filter) is identical to a 
continuous wavelet transform with the change of variables s —►l/s 

2.4 Narrowband Derivation 

For many years wideband processing was considered too computationally expensive. As 
such, an approximation was made so that the processing could be carried out using a Fourier 
Transform instead of a Wavelet transform. In this processing, the time-bandwidth product and 
the velocity are further restrained. Assuming the transmitted signal is of the form of a pure tone 
of center frequency (D c modulated by a complex envelope function w(t), f(t) can be written as: 

f{t) = w(t)e ia,ct (8) 

The time delayed and time scaled received signal, g(t) then becomes 

g(t)=w(s[t-T])e imA ‘- T] (9) 

Defining.? = 1- 8 and 8 - 2 v/(c+ v) , g(t) can be written as: 

g{t ) = wi[l- S][t- T])e im ‘ l, - T] e^ sl ‘- T] 


( 10 ) 
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Defining the Doppler shift along the line of sight to be 0) d = (O c 8 , and assuming v/c«l, 
yields: 

g(t) * w(t - r)e ia>At - T] e i0,At - T] (11) 

or, 

g(t ) » /(/ - r)e ,mAt ~ T] (12) 

This is the narrowband model for signal reflection. The signal has been delayed by x and 
shifted in frequency by an amount co d . It can be shown [29] that this constrains the time- 
bandwidth product such that: 

TB « (13) 

2v 

Again, in order to detect the range and velocity of an object, the received signal g(t) must 
be correlated with hypothesized scaled and time shifted versions of the transmitted signal. The 
m a ximum of the magnitude of this narrowband correlation will correspond to the velocity and 
displacement of the object. 

2.5 Wideband and Narrowband Differences: an Example 

Consider a linear frequency modulated signal with center frequency 50Hz, and bandwidth 
of 80Hz (frequency changes from 10Hz to 90Hz). If this signal were frequency shifted by 100Hz, 
its time duration would be the same, its center frequency would be 150Hz, and its bandwidth 
would remain 80Hz (frequency sweeps from 110Hz to 190Hz). If, instead, this signal were 
scaled by a factor of three, its time duration would be one third of the original, its center 
frequency would still be 150Hz, however its bandwidth would be 240Hz (sweep from 30Hz to 
270). This is shown in figure 1. This shows that scaling a signal yields a different time series than 
when applying a Doppler shift to that same signal. Thus, it is important to process wideband 
signals with time-scale methods rather than with narrowband techniques. 


11 



Figure 1 A comparison of frequency shifting and scaling. A)LFM 
with center frequency 50Hz and bandwidth of 80Hz. B)Frequency 
shifted by 100Hz. C)Scaled by a factor of 3. 


2.6 Deriving the Likelihood Ratio 

Now that it has been shown how the scattering process affects the transmit signal, a 
detector to estimate the scattering processes can be derived. Basic detection and estimation 
theory has been thoroughly laid out in books such as Van Trees [30], [31] and McDonough and 
Whalen [32]. This derivation follows the one in Van Trees [30] chapter 2. In an active system, a 
known signal, f(t), is transmitted into a channel and reflected by scatterers and received as, r u (?) : 

r u (0 = y u ( 0+«*(0 (I 4 ) 

where n u (t) is unwhitened noise, and y u (t) is that part of the received signal caused by the 
reflection and reverberation of the transmit signal. The subscript u denotes the unwhitened input 
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signal. The Likelihood ratio is: 


A ( r ) = 


Pr\H a ( r l-^o) 


where H, denotes the hypothesis that a scatterer is present at a particular scale and delay, and H 0 
denotes the hypothesis that there is no scatterer present. The expression p r ^ H (r\H i ) denotes 
the probability that the observed vector r occurred given hypothesis H, is true (i=0,l). When the 
components of the vector r are jointly Gaussian the likelihood ratio can be written 


A(r) = 


*0 

exp 

(r H - rriy)R { \r- m x ) 

R, 

exp 



(16) 


where and R x are the covariance matrices for the received signal under H 0 and H[ 
respectively, and m 0 and m x are the mean vectors of the received signal, H represents the 
Hermitian transpose (complex conjugate transpose), and rj is the threshold. Taking the log, 
assuming the means are zero( m 0 = n\ = 0 ), and moving constant coefficients to the right hand 
side, the log likelihood ratio is 

l(r) = log A(r) = r H (R~ l - R; 1 )/ y (17) 


where y is 77 divided by all the constant coefficients of equation 16. For the signal defined in 

equation 14, the covariances are: 

Rr. = R„ „ and R, = R r _ = R v „ + R„ „ 

v 1 r u r u y„yu n tt n u 


Here R,, „ is the unwhitened noise covariance, R„ „ is the total echo covariance, and R r , is the 

'hr‘u y u y u r u r u 

received signal covariance. Substituting these into equation 17 yields: 


l =r (Kn, 



Rearranging yields 


/ = r H R:\ 
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Usin 8 R ~ r u r u = R y.y. + 

l= rHR rl( R y,^ R n^~K u OK\/ 


which can be simplified to: 


l=r R r,r. R y.y. R n.n. r 


Identifying that the estimate of the target signal y is: 

y = K. R y.y. r 


( 21 ) 


the likelihood ratio can be written: 


i = y 


D-l 


Knf 


This detector is referred to as the estimator-correlator, because it first makes an estimate 
of the target signal ( y ) and then correlates it with the whitened received signal ( R~^,J )• This 

was first brought to light by Price [33] for Gaussian signals in Gaussian noise. Work has been 
done to extend this to other models of continuous signals by Esposito [34] and Kailath [35], [36], 
and discrete signals by Schwartz [37], [38]. Middleton [39] has analyzed simultaneous detection 
and estimation 

2.7 The Whitening Process 

To derive a detector for a signal in colored noise, a detector is proposed (without 
explanation), and then shown to have suitable detector qualities. This derivation follows the one 
in Van Trees 30 section 4.3.1. Consider a received signal, r u (t), which consists of a desired signal 
y u (t), and colored noise, njt) 

y u (t)+n u (t) 

The filter h w (t,u ) which is known to whiten n^t), is now applied to r u (t) 

r(t)= \h w {t,u)r u (u)du 
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The subscript u is dropped from whitened variables. Expanding yields: 

r(t)= \h w (t,u)y u (u)du+ \h w (t,u)n u (u)du 

The integrals are now renamed 

r(t) = y(t)+n(t) 

Thus the problem has been reduced to detecting a known signal in white noise. The sufficient 
statistic A for a known signal in white noise is: 

lnA[r(0] = \r(t)y(t)dt-^\y 2 (t)dt 

The whitened known signal y(t) can be pre-computed before the likelihood ratio is computed 

so that only the first term on the right hand side needs to be computed in real time. However, in 
this formulation the received signal must be whitened in real time for the likelihood ratio to be 
computed. This is known as the Whitening realization. Another realization can be derived which 
does not require whitening the received signal in real time. Substituting in the expression for the 
whitened variables: 

In A [r(/)] = ldt\duh w (t, u)r u ( u) J dv h w ( t,v)y u (v) - E s 

where E s is the energy in the transmit signal. Changing the order of integration (which can be 
done because the signals are of finite duration) yields: 

In A [/*„ (0] = \dur u (u)^dvy u (v)\ h w (t, u)h w (t, v)dt 
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Defining: 

Q(u,v)= jh w (t,u)h w (t,v)dt 

Substituting this into the log-likelihood yields: 

lnA[r(0] = \dur u {u)\y u (y)Q{u,v)dv 

This can be further simplified by defining 

g(u)= f y u (v)Q(u,v)dv 

which yields 

In A [>(*)] = \r u {u)g(u)du 

This is called the correlator realization, or the mismatched filter for colored noise. It is used 
throughout this work. Because g(u) can be pre-computed, this realization has the same real-time 
computational complexity as standard detection in white noise. 

2.8 Calculating the Eigenvalues of the Noise Covariance 

The inputs to the likelihood ratio (equation 14) are the received signal covariance and the 
noise covariance. It is now shown that the noise covariance can be calculated from the scattering 
function and the transmit signal autoambiguity function. The following development comes from 
personal communication with L. H. Sibul. For a good review of eigenvalue decomposition see 
[45]. The eigenvalues X k and eigenvectors u k (t) of the unwhitened noise covariance are: 

•WO- 
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The noise covariance is composed of white noise and a component due to unwanted scatterers 
and reverberation in the channel: 

n u (t ) = n w (t) + reverb(t) 

where n w (t) jointly Gaussian white noise, and reverb(t ) is the signal reflected by the 
reverberant scatterers. Therefore can be decomposed as: 

ry+Xreverbit’T) 

Where (N 0 / 2 )S(t - r) denotes the white noise component, and Ry^^it, r) is the covariance 
of the reverberant signal. The eigenvalues now become: 

>W0 = f 

The eigenvalues are now defined in terms of a white and reverberant component: 

AT 

1 — __JL_l 2 rev 

A k ” 2 +/i * 

The reverberant eigenvalues are then: 

\ rev u k (t)= \R rev (t,T)u k (T)dT (40) 

The reverberant time signal can be represented in terms of the reverberant scattering function, 
from Sibul, et. al. [1] as: 

reverb(t) = \w rev {x)[U{x)f{t)]d^{x) (41) 

where w rev (x) is the narrowband or wideband reverberant spreading function, dfj.(x) is the left 

invariant Haar measure on the appropriate (narrowband/Heisenberg or wideband/affine) group G, 
U(x) is a unitary representation of the corresponding group, f(t) is the transmit signal, and the 
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vector x represents either (r, <j>) or (r, s ) where x is time delay, <f> is frequency shift, and s is 

scale. The spreading function w(x) can be estimated by wavelet transform of the reverberant 
signal. Equation 41 is essentially representing reverb(t) in terms of the inverse wavelet transform 
of the spreading function. U(x) is an operator which represents the transformation 
f(t) => f(s(t - r)) . A good review of a group theory approach to the ambiguity function can 

be found in Grove [18] or Chaiyasena [40]. The covariance of the reverberation is: 

R rev (t,x)= E{rev(t)rev\x)} 

Rewriting this is terms of the scattering function yields: 

| \E{w(x)wUx')}[U(x)f(t)}{U(x')f(T)}-d^x)d f 4x') 

Defining the reverberant scattering function as S rev (x,x ') = E {w(x)w *(jc')} yields: 

j \s r Jx,x')[U(x)f(l)][U(x')f(T)]-d / u(x)d M (x’) 

Assuming uncorrelated scattering (WSSUS), i.e., S(x,x') = S(x,x)S(x- x') , the 
reverberation covariance becomes: 

\s m (x,x)[U(x)f(t)]lU(x)f(T)]‘dn(x) 

Substituting this into the eigenvalue equation 40 yields: 

WO - f I S r Jx, x)[U(x)f(t)][U(x)f(j)U t ( T)drdMx) 

Taking the wavelet transform of both sides with respect to f(t) yields: 

xr\iu(x')m]-u k (t)dt= 

11 \S r „(x,x)[U(x)fmU(x')f«)]-[U(x)f(T)U k (t)dTdtd»(x) 
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Using the definition of the reproducing kernel for the space spanned by (g, U(x)f), we 
identify: 

K(x,x')= {U(x)f,U(x')f) H 

where(• ,*) H denotes the inner product on the Hilbert space H. K contains all the information 
about the transmit signal and its autoambiguity function. The eigenvalue equation becomes 
X"W, «»(*') = jS r Jx,x)K(x,x'W / u t (x)dp(x) 

where W f denotes the wavelet transform with respect to f(t). The above can be discretized by 
only considering the points ( Xj ) where scatterers (highlights) exist. 

X?W f u k (x t ) = X K{x l ,x j )S rev {x p x J )W f u k (x j ) (50) 

j 

where Xj represents the scale and delay of each highlight and S( Xj , Xj ) represents the strength 

of each highlight. In section 3.6 it is shown that this discretization may dramatically reduce the 
computation time that might have otherwise been needed. Using y/ k = W f u k ( Xj ) the 
reverberation covariance eigenvalues can be written as 

XTft = KSrJ'k (51) 

These eigenvalues and eigenfunctions will be used to whiten the received signal covariance and 
target scattering signal covariance. Similarly, the eigenvalues for R y y are: 

K9k = KS^y? k 

where S laiiel denotes the target scattering function . 

The scattering function of the reverberation and of the target are physical models of 
targets of interest and the environments in which they might be found. The likelihood function 
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from the scattering function, the following holds: 

R" = U y A y Uy + 7 


By noting 


U IU H = U U H = I 

y y y y 


The identity matrix can be brought inside 

R„ = U y [\ y + I]U? 


The likelihood ratio derived in equation 14 was 


l = r H 


Substituting in expressions for the whitened variables yields: 

/= r H [I][U y RyU?W y {R y + iy l U?]r 


which can be simplified to: 

/ = r H U y \ y (\ y + l)~ l Uyr 


Defining 

Z = UyT 


the likelihood becomes: 

i=z H A y (\ y +ir l z 


The covariance of z is now computed for later use: 
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Expanding^ yields: 

K = UfiRy, + l]U y 

or more simply: 

R zz = A y + / 

We have now come quite a ways from our original expression for the likelihood function. We will 
now condense together only the terms needed to calculate the likelihood ratio in practice. First, 
the likelihood function itself is: 

/ = z H \ y (h y + iy'z 


where z is the whitened received signal: 


z = U y r 


U y and A y are the eigenfunctions and eigenvalues of the whitened target scattering signal 
covariance which is: 


R vv = R- V2 U"R vv U n \ n 

yy n u n u y u y u n u 


- 1/2 


U n and A n are the eigenfunctions and eigenvalues of the unwhitened noise covariance which are 
derived from the reverberation scattering function and the transmit signal autoambiguity function: 

= gSLr* 

Likewise the unwhitened target scattering signal covariance R is computed from the target 
scattering function and the transmit signal autoambiguity function: 


Thus the inputs to the likelihood function consist of the received signal, the transmit signal auto- 
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Figure 2 Ingredients of the Likelihood Ratio 


ambiguity function, the target scattering function, and the reverberation scattering function. This 
is shown graphically in figure 2 . 

2.9 Calculating the Receiver Operating Characteristics 

To evaluate the performance of the likelihood detector that has been derived in the 
previous section, it is necessary to find the Receiver Operating Characteristics. The following 
derivation of the ROC comes from personal communication with L. H. Sibul. First, the 
characteristic function of the probability density of the likelihood function is found. The 
probability density will be evaluated by taking a Heaviside expansion of the characteristic 
function. Finally, the probabilities of detection and false alarm will be computed from the 
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probability density function. The log likelihood can be written in scalar form as: 


/ ¥ 

/=! 


i* 4 


^+1 


where N is the number samples in the received vector (as explained in section 3.1), z k are the 
individual samples in the whitened received vector z , and A y = diag(A y k ) . Using 

£{|z,| 2 |«»}=land£{|z 1 | 2 |/?,}=^ + l 
Unit variance white Gaussian random variables can be written in terms of z k 

v * 0 = 7 " 311(1 v *. = k = l ’-’ N 

where sub-subscript 0 or 1 indicates that this variable corresponds to Hq or H, respectively. 

Using these, the likelihood ratio can be rewritten such that the magnitude of each term in the sum 
is expressed solely in terms of the eigenvalues 


.P-yfr-o l 

k =1 A k + 1 *=1 


Using i=0,l, and rewriting the complex variable v k = c kj + jd k yields: 


- 1 K 2 + O ITT; and /|ff, = f (c t , 2 + 




Each term in the sum is of the form = A k x 2 where: 


A k = * ~ and A k = X y k 

*° JZ+1 *' * 


We may now use the property that any random variable that is the sum of independent random 
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variables, has a probability density which is the convolution of the probability densities in the sum. 
This can be found in standard texts such as Papoulis [41], or Bendat and Piersol [5]. In other 
words, for general random variables x, y, z with probability density fiinction p x , p y , p z , if. 

z{k) = x(k) + y{k) 


then 

p(z)= \p x {x)p y (z-x)dx 

or 

where <f>{co) is the Fourier transform of a probability density, and is called the characteristic 
function: 

(/>{(») = E[e Jm ]= fp(x)e ja *dx 

The Characteristic function of one term of the likelihood ratio is: 

(f> y = '* 1,2 dx or </> y = (1 - j2o)A kj )“ 1/2 

Using the property that functions convolved in the probability domain are multiplied in the 
characteristic domain yields: 

^ = n (>- ft<aA t y m (l- j2mA K Y m 

k=\ 


or 

4-n (i-^r 1 

*=1 
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By using the Heaviside Partial Fraction Expansion [42], this can be written as 


N 

A = X ~ 

*=i 1 


K„ 


j2aiA ki 


( 82 ) 


where 




J =1 


1 

(1 - j2oA J ) 


co=l/U2A kj ) 


Subsequently: 


j v A - 


This integral is often listed in Laplace transform tables, such as Spiegel [43], so the following 
substitutions are made: s = jet ), ds - jdco ,dd) = ds/ j 


A Ire", 


Letting JL denote the Laplace transform yields: 


P, 


*=l 


1 - 2sA 


or 


P,,„,(/!#,)= I K k l' 


k= 1 


-2A k js-\!2A k j) 


Evaluating the inverse Laplace transform yields: 


P, <H <jm=± 4 r e ~" 1A ‘' 

*=1 


( 88 ) 
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The probability of detection is defined as: 


P d = |p w (/ \H,)dl 

Y 


Where y is the threshold from equation 17. Substituting equation 88 into the above yields: 


= ry 


s, k l 


-1/2 A t 


^ M ~ 2A k, 


dl 


or 


1 


K k x r -U2A^ 

-2 A k v 


Integrating over / yields: 

A K k[ e ,nA - 

*=1 ~ 2A k t ~^^ 2A k x 

or 

^ = Z K t/ rl2A ‘' (,3) 

yfe=l ' 

Similarly, the probability of false alarm is: 

P f = t W 

Ar=l 

A ROC curve is the plot of probability of false alarm versus the probability of detection. 
Equations 93 and 94 are what are used to compute all of the ROC curves in this thesis. Because 
of the form of equations 93 and 94 two scatterers that have the same amount of energy in the 
whitened eigenvalues will have different ROC curves if the energy is distributed differently 
between the eigenvalues. This is why a scalar figure of merit can not be used to characterize the 
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WTD-EC. In short, the performance can not be characterized by a scalar because the scatterer 
can not be characterized by a scalar. 

One of the original contributions of this thesis is that these curves are not using equal 
eigenvalues as in Roan and Sibul [4]. 
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Chapter 3 

IMPLEMENTING THE WAVELET TRANSFORM DOMAIN 
ESTIMATOR CORRELATOR 
3.1 Implementing the WTD-EC 

One decision to be made in implementing the WTD-EC is deciding how large the matrices 
will be. It is easiest to start at the end. When the detector is actually implemented it will have to 
compute: 

/ = r H R~lRR- x r 

nn yy rr 


In order to include all of the energy reflected from the target, the received vector r has to be at 
least as long as the transmit signal convolved with the scattering function. We will call this length 
N. Now let’s go to the beginning. To find the unwhitened covariance eigenvalues, we must 
implement the equation: 

4</,(*,) = Z K(Xi,Xj)S(Xj,Xj)y/ t (Xj ) (96) 


where each Xj is a position in scale and delay space Xj = ( Sj , r y ). There is no restriction on 

where these positions should be, or what order they should be in, only that there should be at least 
N of them. Again, in order to capture all of the energy in the target, there should be an x j for 

each highlight. 

Theoretically it is possible for there to be more than N highlights (if there are many 
highlights at many scales). In this case it would be necessary to make the received vector longer. 
Equation 96 is easily implemented as a matrix equation: 

= XSy, (97) 

In our case, this equation is implemented with the Matlab eigO function. Matlab will 
automatically normalize y / t . In practice there is one large eigenvalue for each highlight, and the 



30 


eigenvector for that eigenvalue looks like the autoambiguity function centered on the highlight 
position. Note that even though x } does not enter into equation 97, one must still keep track of 

it in order to take the inverse wavelet transform of y/ i . 

W , = Wj'yr^Xj) 

or 

«/ = Z W{Xj)f{t)]y/XXj) 


so u f is composed of scaled and delayed versions of the transmit signal, weighted by y/ i (Xj ). 

Again, this necessitates u being at least N long. Also, the vectors u t are not unit normalized. 
Next, the unwhitened target signal covariance is computed 

R yr =U\ U~ y 


where U v is a matrix whose columns are the eigenvectors u. . Because w, are not unit 
y if • • 

_] JT 

normalized it is not true that U„ = U,, . Also, because its condition number is infinite the 

7U SU 

inverse can not be found in the normal way (such as the Matlab inv() routine). In practice it is 
found that it is nearly true that U V U*? - cl where k is a positive real number (this is shown in 

section 3.3). In the code used for this thesis we use t/J 1 = (1 / c)U**, and similar expression 

su v 7 w 

for the noise covariance eigenvectors. 

3.2 Input SNR 

Throughout this work the input SNR will be used to quantify the received signal seen by 
the detector. It is not clear how to apply the traditional definition of the input SNR to a scatterer 
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with more than one highlight. The traditional definition of Input SNR is: 

E 

InputSNR - — 

E„ 


( 101 ) 


where E s is the energy in the transmit signal reflected by the target, and E n is the noise energy 

integrated over the duration of the received signal (this duration is called the processing window). 
The difficulty is that when a target is spread, the length of the received signal will be longer than if 
the target can accurately be modeled as a point scatterer, and will only have the same input SNR 
if the processing window is increased for the spread target, and not for the point target. This is 
not necessarily a fair comparison. 

One can rewrite the energies in equation 101 as an average power times the processing 
window duration. The processing window duration in the numerator and denominator then 
cancel out, and you get a definition of the input SNR which is independent of the processing 
window 


InputSNR = 


A 

P n 


This is the definition preferred by Scharf [44]. This definition has the advantage that, when the 
energy in a target is spread among several highlights, the average power goes down while the 
noise variance remains the same, causing the input SNR to decrease. This reflects the fact that 
the matched-filter output SNR would also decrease due to the splitting loss incurred by having 
several highlights. 

In terms of the scattering function, the signal power is equal to the magnitude of the 
average target scatterin g function highlight divided by the number of sample points in the transmit 
signal. This is because the magnitude of the scattering function at a point is the total energy 
reflected by that highlight. Because the highlight is reflecting the transmit signal, that energy is 
distributed across the duration of the transmit signal. 

In all of the examples in this thesis, the reverberant scattering function is equal to one at 
each sample point. Because each reverberant highlight is reflecting the transmit signal, it is 



32 


contributing noise to all the sample points around it. For this reason the noise power is also one. 
Note that if the target scattering function also had a large number of highlights very close 
together, its signal power would have to be found in this way. 

Combining these definitions, the input SNR is InputSNR = 10 log / (A^s^)] 

where s'^ is the magnitude of the largest target highlight, N s is the number of sample points in 


the transmit signal, and s^ g is the average reverberant scattering magnitude. 101og, 0 is used and 


not 201og 10 because the scattering function is a second order statistic. Note that if you multiply 
the target and reverberant scattering functions by constants the resulting input SNR is: 


InputSNR = 101og 10 


(<Q 

(bsL)K. 


101ogJC/(C^)]+101og 10 [a/6] 


The signal powers corresponding to a target and reverberant scatterer, and transmit signal 
with 50 sample points are shown in figure 3 (only the real part of the complex envelope is shown). 
A single target highlight of magnitude 50, results in a target signal of power one that is 50 sample 
points long. A reverberant scattering function which is one everywhere, produces a reverberant 
signal with power equal to one everywhere. Therefore, both have the same total energy in both 
the scattering functions and time series. 
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It can be shown that one need not rerun the entire calculation if the only change is the 
input SNR. If the scattering function in equation 97 is multiplied by a constant 

4V, f = K(aS‘)tf and = K(bS rev )v7 


the new eigenvalues are 

/t; = aX° and X\ ev = bX™ 

where X° and X t ev0 are the old eigenvalues, that is, the eigenvalues that would result if one had 
not multiplied the scattering functions by a constant. The eigenvectors remain the same. 
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Subsequently, the whitened covariance matrix is: 

- 1/2 T T H l 


R„ = (*A u“ [£/„ (aA „ )U« ]£/„, (AA ) 


- 1/2 


or 


R„ = (db)R°“ 


and therefore 

A y = (a/b)K yo 

Thus if you want to know the ROC for a given scattering function and transmit signal, but for a 
different input SNR, simply multiply the old eigenvalues by the ratio of the constants and generate 
a ROC curve from the resulting eigenvalues. 

3.3 Numerical Problems Computing the Heaviside Expansion 

The probability of detection and false alarm depend directly on computing a Heaviside 
expansion of the characteristic function of the probability density. While the Heaviside Expansion 
is an elegant mathematical procedure on paper, it can be difficult to implement in practice. It 
turns out that there can be numerical problems in computing the probability of false alarm. It will 
be shown that this only happens when the number of target highlights approaches the number of 
digits of precision of the computer doing the calculation. When this occurs, the probability of 
false alarm can be approximated with a chi squared distribution. Details are given in Appendix B. 
All of the examples in this thesis have a number of highlights less than this limit. 

The probability of false alarm depends on K ka which in turn depends on A k . Because 

A ka = / (A k + 1) , when the eigenvalues are much larger than one, all of the A ko ’s are 

approximately equal to 1. This can lead to subtracting two numbers which are close together in a 
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denominator, a situation that introduces numerical noise into the calculation 

K> " = 0 1- ih 


,,k 


(0=l/j2A t(J 


or 


K . n_!- 


For almost all A, the term in the denominator is 1 minus a number close to 1. This can be 
*0 

avoided by writing A kj in terms of the eigenvalues and rearranging. 

% 

A, -~ ^+i 


*,=n 

j* k 


n{ 

l- -—- 


j2X\ UX y ) 


^-wri 

]*k 


\+X y 


\K-x y jj 


(in) 


This alleviates the problem of subtracting two close numbers in the denominator, and this is how 
the Heaviside coefficients are computed for the examples in this thesis. However, equation 111 
reveals a deeper problem. liX y k is of order one (of which many are) then K ko will be of order 

10*, where N is the number of eigenvalues. This is a problem because the K ko ‘s must sum to 1. 
Because ^ (at) is a characteristic function, it must be true that <j)^ (ct) = 0) = 1. 
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Substituting this into equation 82 yields: 


N 


K t 




I 



Since the K, s must sum to one, the computer must be able to represent a difference of 1 
between two numbers of order 10* . Therefore N must be less than the number of digits of 
precision of the computation. This limits the number of eigenvalues to the number of digits of 
precision of the calculation. As will be shown in section 3.3 the number of significant eigenvalues 
is approximately equal to the number of target highlights. Therefore, in order to compute ROC 
curves for a target with more than N highlights an approximation to the Heaviside function must 


be used. 


This problem persists in other ways of writing the Heaviside expansion, such as: 

-U2A k r l 


N 


a =n 


*=i (J2A ki ) ) 


The Heaviside Expansion would then be 



Q 

(' ®-(j2A ki y l ) 


where 


v 


C,=(a>- 0'2^)-')n 


■C/2 A)' 1 




e)=02^r 


(116) 
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and 


N 

4«>)=I 


N 


£,(0)-C/2 4) “ 


-=S-i2^Q, = i 


where 







W + 1 ) 


( 118 ) 


Again the coefficients are of the order 10 N1 . 

It is important to note that computing the probability of detection does not have this 
problem and can be computed exactly. This is because A kj = k\ and thus they are not 
necessarily all nearly equal. This allows us to find an approximation to the false alarm that does 
not necessarily have any information about the eigenvalues in it. The ROC will still contain 
information about the eigenvalues because they are contained in the probability of detection. 

If all of the A t ’s were precisely equal, the corresponding probability function would be a 

*0 

Chi Squared distribution. This method is not used in this thesis, but may be useful in the future, 
and so details are given in Appendix B. 

3.4 Validating the numerical calculation 

One of the contributions of this thesis is the numerical computation of the Wavelet 
Transform Domain Estimator-Correlator from the transmit signal and scattering functions when 
the covariance eigenvalues are not all equal. In this section an analytic special case is derived to 
validate the numerical calculation. Consider the case where the autoambiguity matrix K is the 
identity matrix I. Substituting into equation 51 yields: 

Xy* = S>' and Xyr* = S>" 

because the WSSUS assumption means that S is diagonal (i.e. S = diag{s k ) ), the eigenvalues 
become the scattering strengths, and the eigenvectors become the usual basis functions 


X k = s[ and A" = S k 



where the ¥ are the matrices formed from the corresponding eigenvectors ¥. Recall that, in 
order to find the eigenvalues from the scattering function, the wavelet transform of the covariance 
eigenvalues was taken. Thus the vectors^* are really the wavelet transform of the eigenvectors of 

the covariance matrices. The matrix of inverse wavelet transformed vectors will be called U . 
Note that U does not equal I, but that U, = U„ . Thus this is not an entirely unrealistic case. 

It is not required that the autoambiguity function of the transmit signal be a delta function, just 
that its autoambiguity function has insignificant sidelobes at delays and scales corresponding to 
the distance between highlights (the Xj from equation 50). The unwhitened target signal 

covariance is: 

Kytyi = U„A U U? (119) 


The whitened signal covariance is: 

R = A~ l/2 UfU,A,UfU„A~ 1/2 

yy n n t t t n n 

But since U„ and U, are equal, all of the U terms cancel and 

Rjy = diag{X / A") 

Since is diagonal its eigenvalues are just the diagonal elements 



( 120 ) 


( 121 ) 


( 122 ) 


The eigenvalues generated by the numerical code used in this paper have been compared to the 
above analytical method for many different scattering functions and they agree to within machine 
numerical precision. 
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For the case where K is not the identity matrix, one can check to make sure that the 
matrices involved conform to constraints on eigenvectors and covariance matrices. To check that 
the inverse wavelet transform of u is done correctly, one can use the property that a matrix of 
eigenvectors is Hermitian. As shown before, one starts with the autoambiguity matrix K and the 
scattering function S, and finds: 

X\j/ = KS y/ 


where the X ’s are the eigenvectors of the covariance matrix and the V ’s are the wavelet 

transformed eigenvectors of the covariance matrix (each of which has either the subscript t or n to 
indicate the target signal covariance or noise covariance). The actual covariance eigenvectors are: 

u - WJ l y/(x) 

Because it is the result of this inverse wavelet transform, the vectors u are not unit normalized. U 
is defined to be the matrix of eigenvectors u , and should be Hermitian. That is 

U H = U~ x or UU H = I 

But because the eigenvectors are not normalized it is actually the case that UU H = cl . The 
product O tj = UU H for a typical set of eigenvectors is shown in figure 4. All of the information 

about the magnitude of the received signal is in the scattering function. The normalization of the 
transmit signal doesn’t matter. Therefore, it may be possible to normalize the transmit signal such 

that c = 1. This would save computation time by not having to compute the product UU in 


order to find c. 
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The covariance matrices R v ,R„ „ and R can also be checked. Covariance matrices 

must be positive definite, which in turn means that their eigenvalues must be positive and real 
[45]. The computation can be checked by seeing how close the eigenvalues come to being 
positive and real. 


Table 1 Closeness of eigenvalues 
to being real and positive 



Negative 

Imaginary 

x rev 

6.62E-02 

0.00E+00 

X 

3.26E-06 

0.00E+00 

x y 

1.19E-17 

3.18E-16 
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Table 1 shows how close the various eigenvalues come to being purely positive and real. The 
“Negative” column shows the ratio of the largest negative eigenvalue to the largest eigenvalue. 
The “Imaginary” column shows the ratio of the largest imaginary component of an eigenvalue to 

the largest eigenvalue. The fact that X y has such small negative eigenvalues is at least in part due 

to the fact that all of the negative eigenvalues were removed from X and X ev before X was 

computed. The largest negative reverberant eigenvalue may not seem clearly negligible in and of 
itself, however the sum of the negative eigenvalues only constitutes 5% of the reverberant 
eigenvalue energy. Negative eigenvalues only constitute 0.01% of the target eigenvalue energy. 
Negative eigenvalues are the result of numerical noise in the eigenvalue decomposition process. 
They must be thrown out because they make the probability of detection and false alarm diverge. 

The code to calculate the WTD-EC is broken roughly into two parts, the first to generate 
the whitened target covariance eigenvalues, and the second to generate the Heaviside expansion 
and ROC curve from those eigenvalues. This section can also be validated. When there is only 
one significant eigenvalue, there is only one term in the characteristic function, so the Heaviside 
coefficient must equal 1 

pr 1 _ 1 


The Heaviside expansion of this is: 



1 

1 - j2mi y 



Eh 

1 - j2oA ti 


Eh 

1 - j2coA { 


meaning K x = 1. The code used to generate the examples in this thesis complies with this 
constraint. 

In [30] pg 108, Van Trees shows that the probability distribution of the likelihood function 
reduces to a Chi Square distribution for the case where all the eigenvalues are the same (equal 
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covariances). His expression for the probability of false alarm is: 


Pf = \~I r 


r 


Vm+ i 




where I r is the incomplete gamma function,is the threshold, and M=N/2 -1 where N is the 

number of eigenvalues. A similar expression is used for P d with a slightly different y"'. The 

code used in this thesis can obviously not compute a ROC curve for eigenvalues that are exactly 
equal, but for a small number of eigenvalues that are close together (in this case within 5%) it 
should yield similar results. A comparison is shown in figure 5. Because only the subroutine that 
takes the eigenvalues as inputs and generates the ROC curves is being validated, the eigenvalues 
used are made up, and do not correspond to any particular transmit signal or scattering functions. 
Note that these curves are not meant to be exactly the same, just to show that no gross coding 
errors were made that might cause the computation to be orders of magnitude wrong. 
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Figure 5 Comparison of ROC computation used in this thesis (exact) and the Chi Square 
approximation 


3.5 The meaning of the whitened target covariance eigenvalues 

As seen in equation 122 the eigenvalues of R ,, v in the analytic case are just the ratio of 
target and reverberant scattering functions at each highlight location. The ROC curve depends 
only on these eigenvalues. Because the ROC is the most valid performance characterization for 
detectors, and the ROC depends only on the eigenvalues, these eigenvalues can be thought of as 
generalized, local, signal to noise ratios. Also, it is easy to see that the number of significant 
eigenvalues is closely related to the number of target highlights. For all the examples given in this 
thesis, the number of significant eigenvalues is equal to the number of significant highlights, 
regardless of the signal or its autoambiguity function shape. It was shown in section 3.2 that if 
the number of significant eigenvalues (and therefore the number of significant target highlights) 
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exceeds the number of digits of precision of the computer, there will be numerical problems 
computing the Heaviside coefficients necessary to compute the ROC curves. 

As we saw in the previous section, when K = I , X y k = s' k / s n k . Recall that 

N 

Therefore, in this case, the likelihood is directly proportional to the sum of the energy in the 
scattering function, regardless of how that energy is distributed. A scatterer with one highlight of 
unit energy should be detected just as well as a scatterer with N h highlights with 1/N h unit energy 
each. Or, a scatterer with N h highlights each of one unit energy should be much easier to detect 
than a scatterer with one highlight of one unit energy. In section 3.6 it is demonstrated that this is 
approximately true even when K ^ / . 

3.6 An observation about the calculation of the covariance eigenvalues 

Earlier, the reverberation covariance eigenvalues were found to be 

^T W f u k( x i) = X K(x i ,x j )S rev (x J ,x j )W f u{x j ) (129) 

j 

where x. represents the scale and delay of each highlight and S(Xj, x j ) represents the strength 
of each highlight. K(x i ,x j ) represents the autoambiguity function at the scale and delay 
corresponding to the distance between scales and delays x ( and Xj . Mathematically speaking 

K(x,x')={u(x)f,U(x')f) H 
K(x,x')= {U(x*x')f,/) H 

where x * x' represents the difference in scale and delay between each highlight. This is the only 
place that the autoambiguity function is used in the calculation of the likelihood ratio. Thus one 
does not need to compute the entire autoambiguity function over the range of interest, but rather 
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one need only compute the points that correspond to the difference in scale and delay between 
each highlight. 

3.7 Properties of the WTD-EC 

The WTD-EC displays many of the advantages that one would expect. For all of the 
examples in this and the next chapter, signals will be referred to by their time-bandwidth product. 
Table 2 shows number of sample points in the transmit signal (N s ), time duration (T), bandwidth 
(B), and time-bandwidth product (TB) of the signals used in this work (for a sampling frequency 
of 4000Hz). 


Table 2 Properties of transmit 
signals 


H 

T (ms) 

B (Hz) 

TB 

201 

50.25 

1250 

63 

201 

50.25 

1093 

55 

201 

50.25 

937 

47 

201 

50.25 

625 

31 

49 

12.25 

1250 

15 

49 

12.25 

1093 

13 

49 

12.25 

937 

11.5 

49 

12.25 

625 

7.5 


Though the WTD-EC and sensitivity analysis developed in this thesis are valid for scatterers 
spread in scale, the computations are only implemented for scatterers spread in delay. Scattering 
functions are sampled in delay with the same sampling period as the signal. Thus the scattering 
function sampling period in distance D s is D s - cT s where c is the speed of sound in the 

medium and T s is the sampling period. For a sampling frequency of 4000Hz the sampling period 

in time T s = 0.25 (ms) . For a speed of sound of c=T500m/s the sampling period in distance 

would be D s - 0.375 m. 
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Figure 6 compares the ROC of one-highlight EC to the WTD-EC for scatterers with N h 
highlights and an SNR of 1/N h compared to the one highlight case. Another way of saying this is 
that the ROC for the one highlight EC for one highlight of unit energy is being compared to the 
WTD-EC for a scatterer with one unit of total scattering strength distributed between N h 


highlights. 



Figure 6 Performance for N h highlights with SNR = 1/N h Signal has TB = 15, highlights are 11 
sample points apart. 


As one might expect, for N h highlights the WTD-EC has very similar performance as the 
one-highlight EC when the highlights have a scattering strength 1/N h that of the one highlight 
case. In a way, it really can be said that it is N h times better for N h highlights. To put this in 
terms of processing gain, one can say that there is a 10Log 10 (N h ) gain in output SNR. 
Additionally, one does not usually consider probabilities of detection smaller than 0.5 useful, and 
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thus the WTD-EC does slightly better in the useful detection range. 

Framing the benefit of the WTD-EC in terms of processing gain leaves out much of the 
picture. Figure 7 compares the ROC for a normal (one highlight) EC with the WTD-EC for four 
highlights. If you consider the one scatterer for the one highlight EC to have a scattering strength 
of 1, the highest curve represents the ROC curve for a scatterer with 4 highlights each with a 
scattering strength of 1, the next curve would be for 4 highlights each with a scattering strength of 
‘/ 2 , and so on. 



Figure 7 Comparing WTD-EC to the one highlight EC. Signal has TB = 15, highlights are 11 
sample points apart, SNRs are Input SNRs. 


In terms of processing gain, for four highlights, the WTD-EC yields a gain of 6 dB. For 
comparison, let’s say that you chose to operate a single highlight detector with a probability of 
false alarm of 10' 3 . This yields a probability of detection of 0.5. Using the WTD-EC would yield 
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a P d = 0.9, which is an 80% improvement in P d . Now let’s say that you chose to operate the 
single highlight detector with aP f = 10' 1 , yielding a P d = 0.8. To achieve a P d = 0.8 the WTD-EC 
could operate at a P f of 1 O' 6 , a six order of magnitude improvement in P f ! This serves to show 
that there are important insights to be found in the ROC curve that are lost in scalar measure such 


as dB processing gain. 



Figure 8 Dependence of ROC on Time-Bandwidth product. There is one highlight, and input 
SNR = -lOdB. 


Figure 8 shows the difference in ROC of the EC between when the noise is correlated, and 
when it is uncorrelated. When the transmit signal is perfectly resolvent (i.e. K = 7) the 
reverberation is uncorrelated. In the section 2.3 it was shown that under these conditions it is true 

that: X y k - s‘ k / s n k . If there is only one highlight, there will be only one non-zero whitened 
eigenvalue, and no need to do a Heaviside expansion. Under these conditions, the probability of 
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detection and false alarm are: 

— -r(i+^') 

P d = e 2Xy and P f = e lx 

which can be further simplified to P f = P d UD , where D is the ratio of the total signal energy to 

the noise covariance. This is the traditional result for the one highlight EC and the matched filter, 
and does not depend on any qualities of the transmit signal except its total energy. When the 
noise is correlated, this level of performance can no longer be attained, and serves as an upper 
bound on the EC performance under correlated noise. As shown, the performance of the EC 
under correlated noise increases with the TB of the signal. 

For every example in this thesis, the reverberant scattering function has a magnitude of 
one (plus a small random number to avoid numerical problems) at each sample point. 

This is quite different than uncorrelated white noise. Having one reverberant highlight at each 
sample point does not lead to having the noise energy evenly distributed between the noise 
eigenvalues the way it would be for white noise. For low TB signals the noise energy is 
concentrated in a few eigenvalues - the energy becomes more and more evenly distributed as the 
TB increases, until the signal is perfectly resolvent, the energy is completely evenly distributed, 
and only then is it really uncorrelated white noise. Thus, all of the examples in this thesis are 
reverberation limited, not noise limited. 
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The point of the WTD-EC is to combine the energy from multiple highlights. Figure 9 


shows that indeed the ROC improves with the number of highlights. 



Figure 9 Dependence of ROC on number of highlights. Signal TB = 15, Input SNR = -lOdB, 
Highlights separated by 10 sample points, and have equal strength. 


When highlights are far apart compared to the width of the mainlobe of the transmit signal, 
the eigenvalues and eigenvectors behave much as if the autoambiguity function was a delta 
function. It would be difficult to predict analytically how the probability of detection changes 
when the highlights are closer than the width of the mainlobe. 
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Figure 10 Dependence of Pd on Separation of Highlights for a transmit signal with TB = 7.5. 
SNR = -3dB. 


Figure 10 shows the probability of detection for a probability of false alarm of 10' 3 for a 
signal with TB = 7.5 and a scatterer with two highlights, and shows how the P d depends on the 
separation between the two highlights. 
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Figure 11 shows the same graph for a signal with a TB = 15. 



Figure 11 Dependence of Pd on Separation of Highlights for a transmit signal with TB = 7.5. 


SNR = -3dB. 

Neither is the a large dependence on the separation, or an obvious correlation with the 
autoambiguity function. 









Chapter 4 

SENSITIVITY OF THE WAVELET TRANSFORM DOMAIN 
ESTIMATOR-CORRELATOR 
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4.1 Deriving the Sensitivity of the Wavelet Transform Domain EC 

Incorporating the estimate of the scattering function into the EC can lead to dramatic 
detector gain. However, if that estimation contains errors, it will not reap the full benefits of 
making the a priori estimate. The critical question is: does the advantage of higher time- 
bandwidth product signals and recombining highlights outweigh the detriments caused by errors in 
the a priori estimate? For this reason, the sensitivity of the EC to errors in the a priori estimate is 
derived. Compelling evidence is given that the WTD-EC does in fact retain it’s advantages even 
when the a priori estimate contains large errors. Preliminary results on this topic have been 
previously published by the author of this thesis [46], [47]. 

An expression for the ROC is now derived that takes into account errors in the assumed 
scattering function. The likelihood function derived before is: 


l=r H U, 


diag 


% 


T\ 


1 + 1 ? 


U H r 

y 


(132) 


For the case where the assumed scattering functions are incorrect, the affected variables are 
designated with the sub/superscript a. The assumed likelihood function is: 

\ 

u H ya r (133) 

J 

Note that we have not made any assumption about how the eigenvalues and eigenvectors are 
found. Thus, the sensitivity analysis derived here is true for any Estimator-Correlator, and not 
just the Wavelet Transform Domain Estimator-Correlator. U ya in equation 133 does not whiten 

r because it is not based on the correct scattering function, r is now re-expressed in terms of a 


l a = r H U 


ya 



r X* 

diag 

l 

\+xr_ 
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whitened vector z a . 


z - A r U r r 


(134) 


It is now shown that z a is white 


R zz = E{z a z aH } 


K = E{A- r y2 U?rr H U r K~ V2 } 


R zz = E{A- r y2 U r H R rr U r A- uz } 


- 1/2 ' 


R zz = {A- r m U?U r A r U?U r A~r} 


Htt A -1/21 


«„ = £{A; 1 , 2 A r A; l,2 }=i 


Rewriting r in terms of z and defining 


D a = diag 




(140) 


the likelihood ratio can be found in terms ofz a 


r = (z-AfufjUj.D.u" (u r Afr) 


(141) 


The terms in the middle can be collected and diagonalized. Defining new variables A ar , and U a 

r H _ a1/2tt^tt t\ tt^ tt a 1/2 


U A U = A U U D U U A 

ax ax ax r r ya a ya r r 


(142) 


l a = z a U H A U z L 

ax ax ax 


( 143 ) 
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This can now be rewritten in terms of a new, unit variance, jointly Gaussian vector x a - z a U^ 

/• = x -trkf (144) 

k 

Since x a is jointly Gaussian one can derive a ROC curve for this likelihood ratio. But first x a 
must be normalized. The expected value of x a is: 

x* = U”z a = U^Ufr 

E{\x°\ 2 }=E{x aH x a } 

E{\x°\ 2 } = E{r H U r K-™UJ/ZA-™U?r} 

E{\x a ?}= E{r H U r kT?U?r} 

Using the definition from section 2.8 

z = Ufr 

and the fact that 

£{|z t | 2 |/f 0 }=l and£{|z t | 2 |/f,}=l+2 t 

yields 

£{KTl^o> = 7r“d£{K| 2 |ff,} = l 

Rewriting x a k in terms of unit normalized values yields: 

< 0 = x a k Xk and wj = x a k 
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and l a \H x = 


I 




a |2 


AT 


The w k terms are unit variance, jointly Gaussian random variables (not to be confused with the 

spreading function w(x)). Thus the likelihood depends only on the assumed eigenvalues AT and 

the actual eigenvalues . This expression includes information about both the correct scattering 

function and the assumed scattering function. Thus one can compare ROC curves for correct 
scattering functions and scattering functions with errors. 

4.2 Estimating the Input SNR 

The scattering function estimate is not just an estimate of where the highlights are, it also 
implies an estimate of what the overall input SNR is. The effect of guessing the wrong input SNR 

can be derived analytically. In section 3.3 it was derived that ifS' 0 = aS‘ then AT = aA y k , and 

that the eigenvectors will be the same for both the correct and assumed covariance matrices. 
Therefore equation 141 becomes 

l a = z aH A y2 D a A y2 z a 


and 


AT = AT 


J±iL 

l+AT 


Making the substitution A k a = aA y k yields: 


AT = a% 


l+aA y k 
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or 



K +1 

Xl + Ua 


If a=l we get the sensible result that Xf = X y k and the performance is exactly the same as if the 
sensitivity analysis had not been done. If a<l the denominator is greater than the numerator and 
X* < X y k . If a>l then the denominator is smaller than the numerator and Xf is slightly larger 

than X y k . Asa approaches infinity the assumed eigenvalues converge to 

LimXf = X y +l 

co 

One should not jump to the conclusion that performance must therefore be better. The X a k in the 
above expressions is really X a k \ H { . X k * | H 0 would be: 

X ya kX y X y 
* 1 °“ l+X ya " 1 +aX y ~ X y + l la 


There is no obvious way to analytically derive the effect on the ROC curve so we will simply 
compute it. Figure 12 shows the effect of overestimating the input SNR. 
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Figure 12 Effect of overestimating the input SNR. Scatterer has two highlights separated by 25 
sample points, signal has TB = 15. 


One can clearly see that the ROC curve is indistinguishably different when one overestimates the 
input SNR. This result was found by multiplying the target scattering function by a constant, no 
change was made to the reverberation scattering function. In the limit as one multiplies the target 
scattering function by as large a constant as numerically feasible, one is essentially estimating that 
there is no reverberation. Thus, in practice, it is not necessary to make any estimate of the 
reverberation scattering at all. This is advantageous in two ways. First, it requires less a priori 
knowledge of the scattering environment. Second, the bulk of the computation time is spent in 
computing eigenvalue decompositions, omitting the reverberation scattering function requires one 
less eigenvalue decomposition. Again, this result is true of any EC, not just the WTD-EC. 
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4.3 Validating the Sensitivity of the WTD-EC 

The simplest way to test the sensitivity calculation is to set the assumed scattering function 
exactly equal to the correct scattering function. The resulting ROC curve should be the same as 
the ROC curve for the correct scattering function alone. This has been done for many different 
scattering functions. 

One can also develop an analytic case similar to the one used to validate the correct 
likelihood calculation. Consider the case where U = U ya - U r . When this is true the 

expression in equation 142 becomes: 

= A. l ' r 2 U H UD a U H UA ll r 2 


or 


U A U H =A m D A 

^ ax lx ax ax r 


1/2 

r 


nax _ nya 


1+Aj 

1 + A? 


Substituting into the likelihood functions under hypothesis H 0 yields: 



i±jC 

rJ 


1 

1+^ 


or 



K 

i 


Notice that this is the same as the expression for the likelihood under H 0 for the correct scatterers, 
except the eigenvalues are replaced with the assumed eigenvalues. 



60 


Substituting the special case eigenvalues into the likelihood functions under hypothesis H, 

yields: 



l±A 

i+Af 


( 161 ) 


This is significantly different than the likelihood under H, for the correct scatterers which was 

k 


As stated earlier, if the assumed scatterer is the same as the correct scatterer then the assumed and 
correct likelihood functions should be the same. If the assumed and correct scatterers are the 

same then the eigenvalues will also be the same X y ° = A y k . In which case the 1 + A, terms in 

equation 161 cancel, and the assumed and correct likelihood functions are indeed equal. Another 
way to say this is that the likelihood is maximized when the guessed scattering function is the 
correct scattering function. 

To test the case where the correct and assumed scattering functions are not equal, one 
needs to create a situation where U = U = U r . Just as in the case for validating the correct 


likelihood function, this can be done by setting K = I . In this case the assumed eigenvalues are 
just the ratio of the assumed target highlights and the assumed reverberant highlights. 



Ja 



As before the correct scatterers are: 
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Again, just as was the case for validating the correct likelihood function, for the case 
where U * U r one can still check the results by checking that eigenvector matrices are 

Hermitian and eigenvectors are real and positive. 

4.4 Results 

The purpose of the sensitivity analysis is to have some method of measuring what types of 
signals are best for what applications. What sort of signal should one use when you know the 
scattering function precisely, or when discriminating between very similar scatterers? What type 
of signal should one use when the scattering function is not well known? If there are sidelobes 
that correspond to the distance between highlights, is this a detriment, an enhancement, or 
irrelevant? Before doing a rigorous analysis, one might reasonably expect that the sensitivity to 
error is related to the autoambiguity function of the signal. That, as the guessed highlight moved 
outside of the mainlobe of the autoambiguity function centered on the correct highlight, 
probability of detection decreases. This has the unfortunate effect of reducing some of the gains 
made by having a high time-bandwidth product signal. As will be shown in this chapter, this is 
true to an extent, but luckily not so much so that errors in the scattering function outweigh the 
benefit of implementing the WTD-EC, or of using higher time-bandwidth product signals. 
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In figures 13 through 15 a simple example case is given to introduce the concept of a 
sensitivity graph, and to show it’s connection to the assumed scattering functions. Figure 13 
shows a set of scattering functions which will be used to generate ROC curves. The first 


scattering function is considered to be the correct scattering function. In successive scattering 
functions the guessed second highlight is further and further away from the correct hi ghl ight 
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Figure 13 Scattering functions used to generate the next two graphs 
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Figure 14 shows the ROC curves corresponding to these scattering functions. The first 
curve corresponds to when the first scattering function is both the guessed scattering function and 


the correct scattering function, the next curve represents the second scattering function being 
guessed and the first is still the correct, the next curve is the third being guessed and the first 
correct, and so on. 



Figure 14 ROC curves corresponding to the scattering functions in figure 13 

Graphs like figure 14 can be very muddled when the sensitivity does not depend monotonicly on 
the error (which is generally the case). To overcome this we will introduce a new kind of graph, 
which will be called a sensitivity graph. A sensitivity graph is generated by taking a graph such as 
figure 14 and slicing it vertically at a particular false alarm. Therefore a sensitivity graph is a plot 
of probability of detection versus error for a particular false alarm. For all of the sensitivity 
graphs in this thesis, the probability of false alarm is 10' 3 . 
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Figure 15 shows the sensitivity graph corresponding to figure 14 



0.55 ! - L ---- 1 -'--—*— -1— -j_ _ 

0 0.5 1 1.5 2 2.5 3 3.5 4 

distance of 2nd hypothesized highlight from 2nd correct (in sample points) 

Figure 15 Sensitivity to error in the position of the second of two highlights for a signal of 
TB = 15 
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The sensitivity generally follows the autoambiguity function of the transmit signal. The wider the 
autoambiguity function, the less sensitive the detector to errors in the a priori scattering estimate. 



Figure 16 Sensitivity to the error in position of the second of two highlights for a transmit signal 
with TB = 7.5, input SNR of -3dB, two highlights separated by 50 points. The autoambiguity 
function is not plotted on the same vertical scale. 
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Figure 17 Sensitivity to the error in position of the second of two highlights for a transmit signal 
with TB = 15, input SNR of -3dB, two highlights separated by 50 points. The autoambiguity 
function is not plotted on the same vertical scale. 


There are many things to note about these two graphs. First of all, the y axis does not extend 
from zero to one on either graph. Thus, even though the lower time-bandwidth product signal is 
less sensitive to error, the higher TB signal has a higher probability of detection everywhere. The 
dashed red line indicates the probability of detection when one guesses only one highlight, and 
uses the same transmit signal. Thus the one highlight detector sees only one of the two highlights, 
and therefore half the total scattered energy that the WTD-EC sees. One might guess, before 
having done this calculation, that the probability of detection for two highlights would converge 
to the performance for one highlight, when one of the guessed highlights is completely outside the 
mainlobe of the autoambiguity function. For this reason, the autoambiguity has been plotted on a 
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vertical scale such that it’s maximum coincides with the probability of detection when both 
highlights are guessed correctly, and it’s minimum coincides with the probability of detection 
when only one highlight is guessed. Thus figures 16 and 17 show that the probability of detection 
is even better than one might have expected. The probability of detection when guessing two 
highlights is always better than guessing only one, even when one is completely wrong about 
where the second highlight is. 

Even though figures 16 and 17 may give us hope, it is not generally true that the 
sensitivity curve of a higher TB signal is everywhere better than the sensitivity curve for a lower 
TB signal. In figure 18 one can see that a signal with a time-bandwidth product of 13 does better 
than a signal with a TB of 15 for errors between 5 and 15 sample points, and is of comparable 
probability of detection for errors less than that. That, however, is not a very large difference in 
TB, and for significantly smaller time-bandwidths the probability of detection is uniformly worse. 
Thus there is no hard and fast rule saying that higher time-bandwidth product is always better. If 
you are trying to discern between very similar scattering functions, then you want the signal to be 
sensitive to errors, and it is probably best to pick the largest possible time-bandwidth product. If 
you want the detector to be robust to errors, then you may want to pick a slightly lower time- 
bandwidth product than the maximum allowed by the available equipment. There is a limit, 
however, to trading TB for robustness to error. One would not choose a time-bandwidth product 
of 11.5 for this scatterer because one could always do better by using a signal with TB of 15. 
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Figure 18 Sensitivity to error in position of the second of two highlights for transmit signals with 
various time-bandwidth products, input SNR of -3dB, two highlights separated by 50 points 


While the oscillations in the sensitivity of the TB = 11.5 detector do, in fact, correspond in 
position to autoambiguity function sidelobes, it is not clear why they are so large. Figure 19 
shows that the TB = 15 signal retains its advantage over signals with TB = 11.5 or lower, even 
when at very low input SNR. 
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Figure 19 Sensitivity to error in position of the second of two highlights for transmit signals with 
various time-bandwidth products, input SNR of-13dB, two highlights separated by 50 points 

Figure 20 shows more directly that the sensitivity does not change dramatically with the input 
SNR 
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Figure 20 Sensitivity to the error in position of the second of two highlights for various input 
SNRs, a transmit signal with TB = 15, two highlights separated by 50 points. 

It will now be demonstrated what happens when one is wrong both about the position and 
number of highlights. Figure 21 shows that probability of detection can only decrease if you 
hypothesize that there are three highlights when there are only two. 
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Figure 21 Sensitivity when hypothesizing that there are three highlights when there are actually 
two. Transmit signal has TB = 15, input SNR of -3dB. The two correct highlights are 50 sample 
points apart. The autoambiguity function is not plotted on the same vertical scale. 


Figure 22 shows the sensitivity when there are three correct highlights, and there are only 
assumed to be two. The three correct highlights are separated by 12 sample points, the first 
hypothesized highlight is correct, the horizontal axis represents the distance of the second 
hypothesized highlight from the first. 
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Figure 22 Sensitivity to the position of the second of two hypothesized highlights when there are 
three correct highlights. Signal has TB = 15, input SNR = -lOdB, correct highlights are separated 
by 12 sample points. 




In figures 21 and 22 the probability of detection is never less than probability of detection when 
only one highlight is hypothesized. In all of the examples in this thesis, when the correct scatterer 
has more than one highlight, it is always better to hypothesize more than one highlight, no matter 
how wrong you are about how many or where they are, than to guess only one highlight. 

Roan [48] has suggested that a boot strap procedure could be used to find the scattering 
function that maximizes the probability of detection for a given observation vector. First, the 
maximum on the scale-delay map is considered to be a highlight. Next one hypothesizes a second 
highlight at each position further and further away from the first highlight. When a maximum in 
the likelihood ratio is found, a second highlight is added to the estimate. More highlights are 
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added to the estimate in a similar fashion. Figure 22 Shows there will indeed be maxima in the 
sensitivity when you get more than one (but not all) of the highlights estimated correctly. Thus, 
this bootstrap procedure has merit. 

There is more to get wrong than the position of the highlights, one can also estimate the 
relative magnitudes of the highlights incorrectly. Figure 23 shows the sensitivity to estimating the 


magnitude of the second highlight incorrectly, for a scatterer with two highlights. Highlights are 
separated by 25 samples points, the signal has TB = 15. 



= 15, highlights are separated by 25 sample points 

This is an encouraging graph. It says that it does not matter if you overestimate, but still gives the 
sensible result that guessing too small is equivalent to guessing that there is only one highlight. 
Combined with the result of section 4.2, it can be concluded that one only need estimate the 
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position of target highlights. 

Thus, the a priori knowledge needed by the EC is much less than is implied by the WTD- 
EC without doing a sensitivity analysis. Even though the mathematical form of the WTD-EC 
includes a detailed target and reverberant scattering function estimate, one need only estimate the 
position of target highlights. The relative magnitude of target highlights, and reverberation 
highlights need not be estimated at all. 

Roan’s thesis [17] tackles the problem of how to track a rotating distributed object. The 
sensitivity analysis can be used to find the sensitivity to error in the hypothesized rotation angle. 
Figure 24 shows how the highlights of a rotating object with one highlight at each end are 
projected onto the line of site of the transmit signal. 



Figure 24 How highlights of a rotating objected are projected onto the line of site of the transmit 
signal. 
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When the rotation angle estimate is too small, the hypothesized highlights are closer together than 
the correct highlights. When the rotation angle estimate is too large, the hypothesized highlights 
are farther apart than the correct highlights. 

The sensitivity to error in hypothesized rotation angle is now computed for an object with 
six highlights. Figures 25 and 26 show the correct and hypothesized scatterers used for this 
calculation. Note that the vertical scale of figure 25 is larger because at an angle of 90° the 
highlights combine to form a single highlight with magnitude six. 



Figure 25 Correct scattering function, and hypothesized scattering 
functions with estimated angle too small. 












Figure 26 Correct scattering function, and hypothesized scattering 
functions with estimated angle too large. 


Figures 27 through 30 show the sensitivity to error in the hypothesized rotation angle 
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Figure 27 Sensitivity to hypothesized rotation angle. Input SNR = -lOdB, 
correct angle = 30°. 
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Figure 28 Sensitivity to hypothesized rotation angle. Input SNR - -13dB, 
correct angle = 30°. 
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Figure 29 Sensitivity to hypothesized rotation angle. Input SNR = -16dB, 
correct angle = 30°. 
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Figure 30 Sensitivity to hypothesized rotation angle. Input SNR = -20dB, 
correct angle = 30°. 
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Plotted in each figure is the probability of detection of a one highlight detector of TB = 2. This 
time-bandwidth product was chosen because its autoambiguity function is wide enough to cover 
all of the highlights. It sees the same total reflected energy as the WTD-EC. Thus the TB = 2 
line represents the best the one highlight EC can do under correlated noise conditions. Therefore, 
even though the WTD-EC is less robust to error in rotation angle at very low input SNR, than at 
higher input SNR, it is still everywhere better than the one highlight EC because the one highlight 
EC is limited to a lower TB in order not to incur splitting loss. Furthermore, there is a clear 
relationship between time-bandwidth product, input SNR, and robustness. 

Figures 27 through 30 were generated by first computing the sensitivity when the 
estimated angle is too small, and then when it is too big. This is why the scattering functions are 
split into to graphs in figures 25 and 26. Small random numbers are added to the scattering 
functions for numerical stability reasons. These small random numbers are different in each run, 
and thus sometimes the sensitivity to the angle being too small, and the sensitivity to it being too 
big, do not quite match in the middle (i.e. the TB = 31 line). 

4.5 Summary 

A sensitivity analysis has been derived which can be applied to any EC. It has been shown 
that it is not necessary to make an estimate of the reverberation scattering function. It has been 
demonstrated that the sensitivity of the EC mimics the autoambiguity function, but that the 
probability of detection does not degrade to the probability of detection of N h -1 highlights when a 
hypothesized highlight is much farther away from its correct position than the width of the 
autoambiguity function. It was shown that one can trade off time-bandwidth product to gain 
robustness to error, but that for a given maximum TB there is a TB below which there is no 
benefit. It was demonstrated that making an estimate of the scattering function is always better 
than not making one (one highlight estimate), regardless of how wrong that estimate is. Finally, it 
was shown that one only needs to estimate the position of target highlights, no estimate of the 
relative magnitude of the highlights need be made. 
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Chapter 5 
CONCLUSIONS 
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5.1 Implementing the WTD-EC 

In Chapter 3 the Wavelet Transform Domain Estimator-Correlator derived by Roan and 
Sibul [4] was implemented. Many implementation and numerical problems were overcome. A 
fast method of inverting the inverse wavelet transformed eigenvector matrices was found. It was 
shown that small negative eigenvalues can arise from numerical errors, and that these must be 
omitted for the probability densities to be finite. A method of constructing and choosing the 
minimum size of the scattering matrix and eigenvectors for a given scattering function and 
transmit signal was presented. 

It was shown that, because of numerical accuracy issues, there is a limit to the number of 
significant eigenvalues that can be used in the Heaviside expansion for the computation of the 
probability densities. The number of significant eigenvalues (which is roughly the same as the 
number of significant highlights) must be smaller than the number of digits of precision of the 
computer. It may be possible to approximate the Heaviside expansion with a Chi Square 
distribution, and some work in that direction is given in Appendix B. 

Various methods of validating the numerical calculation were derived. It was shown that 
if the autoambiguity matrix is the identity function then the eigenvalues are the ratios of the target 
and reverberant scattering strengths. One can verify that the inverse wavelet transformed 

eigenvectors were correctly computed by verifying that UU — cl . Because covariance 
functions are positive definite, and positive definite matrices must have real non-negative 
eigenvalues, one can check the accuracy of the computation by seeing how close the actual 
computed eigenvalues come to meeting this condition. Lastly, the subsection of code that takes 
the eigenvalues as inputs and computes the probability distributions can be checked by giving it N 
eigenvalues that are close to the same value, and comparing the result to a Chi Squared 
distribution with N degrees of freedom. It was shown that the code used to generate the 
examples in this thesis meets all of these conditions. 
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Finally the properties of the WTD-EC were presented. It was shown that for N highlights 
the WTD-EC yields a 10Log 10 (N) dB processing gain. Or, more precisely, the ROC curve for an 
N highlight target where each highlight has a scattering strength of 1/N is nearly identical to the 
ROC curve for one highlight with a scattering strength of 1. This is why one can not use just a 
scalar figure of merit to characterize the performance of the detector. Scatterers with the same 
total scattering energy can result in different ROC curves. The WTD-EC also has the expected 
properties of better performance for higher time-bandwidth product signals, and better 
performance for a larger number of highlights. 

5.2 Sensitivity Analysis of the WTD-EC 

A method was derived for finding the ROC curve of the EC for the case where the a priori 
estimate is not correct This result is not only true for the wavelet transform domain 
implementation of the EC, but for any method of finding the covariance function eigenvalues. 
This analysis has provided powerful evidence that the gains of the WTD-EC derived in Chapter 3 
can actually be achieved in practice and provides a method of designing signals for specific 
applications. 

Methods for validating the sensitivity calculation similar to those used in Chapter 3 were 
applied and met. 

A new performance measure was introduced which we chose to call the sensitivity curve. 
The sensitivity curve is the probability of detection versus some measure of the error in the 
scattering function, at a given false alarm rate. It was shown that the sensitivity curve generally 
follows the autoambiguity function. That is, signals with a high resolution are more sensitive to 
errors in the scattering function estimate. This does not mean, however, that there is no 
advantage to using a high time-bandwidth product signal. Because the overall performance 
decreases as the resolution decreases, there is a limit to how much one can trade off resolution for 
robustness to error. For a signal with a sufficiently lower resolution, its performance will be 
worse no matter how large the error is. If you are trying to discern between a set of very similar 
scatterers, and you do not want robustness to error, it is always best to use the highest time- 
bandwidth product signal possible. 

It was shown that the sensitivity does not depend on the input SNR. 
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Most importantly, it was shown that it is always better to make a scattering function 
estimate than not to, and less a priori information is needed than one would expect. It is not 
necessary to estimate the reverberation scattering function, or the relative magnitude of the 
scattering function highlights. 

It was concluded that the WTD-EC can result in dramatic gains over the one highlight EC, 
and that these gains are not just theoretical. Errors in the scattering function estimate do not 
erase the gains made by using a high time-bandwidth product signal. 

53 Future Work 

Many questions remain that can be answered with the theory and computational code 
developed in this thesis. In Chapter 4 it was shown that if one guesses that there is an additional 
highlight which is close to a correct highlight, the performance will not be significantly degraded. 
Can one increase the robustness to error by doubling up highlights (that is, replace single 
highlights with several highlights packed together) ? Most of the examples in this thesis are for a 
small number of highlights. Some examples with a larger number of highlights has been computed 
by Ricker [19], but this should be demonstrated for a greater variety of examples. Also, does 
anything significant happen when a scatterer has some highlights that are within one mainlobe s 
distance of each other, and some that are not? Our derivation did not assume discrete scatterers, 
but all of our examples have been composed of discrete scatterers. What can be said about 
continuous scattering functions? 

Most of the results in this thesis are empirical, they result from computing the performance 
for many combinations of transmit signals and scattering fiinctions. It has been shown that the 
likelihood function increases in proportion to the energy in the scattering function (regardless of 
how it’s distributed between highlights). In chapter 4 it was demonstrated that the sensitivity is 
related to the autoambiguity function. It would be nice to show this analytically. This would 
require either an exploration of the eigenvalue decomposition of the covariance matrices, or 
expanding those matrices in terms of a predetermined set of basis functions. 

There is much room to improve the computational implementation of the theory. Though 
the WTD-EC and the sensitivity were derived for scatterers spread in both delay and scale, the 
computation was only implemented for scatterers spread in delay. The bulk of the computational 
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time is spent finding the eigenvalues and eigenvectors of the covariance matrices. This is done 
with the Matlab eigO function, which makes no assumptions about the form of the matrix. In this 
thesis though, quite a lot is known about the matrices. For example, in: 

hVk = My* 

it is known a priori that S is diagonal, that K is composed of translated versions of the 
autoambiguity function, and that K is Hermitian. The whitened covariance matrix is composed of 
the eigenvalues and eigenvectors of the scattering and noise covariances. Certainly with all this a 
priori knowledge there is some more efficient way to compute the eigenvalues. 

There is fertile ground for new research by relaxing some of the restrictions imposed in 
this work. In the very beginning it was assumed that the highlights in the scattering function were 
uncorrelated. 

S(x,x') = 5(x,x)^(x- jc') 

This made a matrix representation of the problem easy. However, there is no reason to assume 
that the computation could not be done without this restriction. 

The WTD-EC can be combined with other methods in useful ways. Roan has developed a 
method of using Kalman filtering to track highlights. This could be used to generate the 
scattering function estimate used by the WTD-EC. 

The WTD-EC can be combined with search algorithms in useful ways. It has been shown 
that the likelihood function is maximized when the guessed scattering function estimate is correct. 
Roan’s boot strap algorithm could be implemented to find the scattering function which 
maximizes the likelihood ratio for a particular observation vector. This scattering function could 
be the first estimate used in the Kalman tracking method mentioned above. 

A search algorithm could be used to find the optimum transmit signal given a particular 
scattering function and a bound on the scattering estimate error. Similarly, a search algorithm 
could find an optimum transmit signal to distinguish between similar scatterers. 

The sensitivity analysis can be used to decide how large a library to build to distinguish 
between the same scatterer in different conditions. The same scatterer under different multipath 
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conditions will exhibit quite different highlight structure. One approach could be to build a library 
of scattering functions composed of the target at different positions in range and depth. The 
sensitivity analysis could tell you how finely you would need to sample in range and depth so that 
the performance is not severely degraded when the target is in between members of the library. 
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APPENDIX A 
Signal Design 

In 1959 Stewart and Westerfield published their paper "A Theory of Active Sonar 
Detection" [50] where the results of radar echo detection are applied to sonar. Stewart and 
Westerfield optimize a detection index, however the detection index is not a complete 
performance measure because it can not be used to uniquely identify the corresponding ROC 
curve. One needs the complete ROC curve to measure the performance of the detector because 
the total scattering strength no longer uniquely specifies the ROC curve. As shown in section 3.6, 
scattering functions with different highlight distributions, but the same total scattering energy 
result in different ROC curves. 

One set of functions that has been optimized for different qualities is frequency hop codes. 
Hop codes have several desirable qualities such as energy efficiency and predictable sidelobe 
levels. They also have a large number of parameters which can be adjusted, namely the time that 
each frequency chip occurs in the signal. Hop Codes are typically designed to have a “thumb¬ 
tack” like autoambiguity function, but in the past the ROC has not been used in the criteria of the 
design. Also, the vast majority of frequency hop code work has used narrowband processing. 
Golomb and Taylor’s “Two-Dimensional Synchronization Patterns for Minimum Ambiguity” [50, 
59] provides a detailed review of the literature up to 1982. 

Titlebaum’s “Time-Frequency Hop Signals Part I: Coding Based Upon the Theory of 
Linear Congruences” [51] describes a code set where there is an upper bound on the cross- 
correlation between any two members of the set. In “Time Frequency HOP Signals Part II. 
coding Based Upon Quadratic Congruences” [52], Titlebaum and Sibul define signal efficiency as 
the ratio of average power to peak power and conclude that maximum signal efficiency is attained 
by frequency hop codes. In this frequency hop coding scheme the distance between successive 
chips increases quadratically. Chips wrap around in frequency when they surpass the desired 
bandwidth. Limits on the side lobe level for Quadratic Congruence codes are derived. 

In “Ambiguity Properties of Quadratic Congruential Coding”, by Titlebaum, Marie, and 
Bellegarda [53] it is shown that these codes have a thumbtack like auto-correlation function with 
a uniform upper bound on the plateau. It is shown that their cross ambiguity properties are better 
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than Costas codes (described below). 

Further work along these lines is presented in Bellegarda and Titlebaums’s “Time- 
Frequency Hop Codes Based Upon Quadratic Congruences” [54] and “Frequency Hop Multiple 
Access Codes Based Upon the Theory of Cubic Congruences”, by Marie, and Titlebaum [55]. 

Bellegarda and Titlebaum, in their 1991 paper “The Hit Array: An Analysis formalism” 
[56], show that discretizing the time frequency representation of a signal provides a convenient 
method for analyzing them. They show that reflections and rotations of hop codes will produce 
the same reflection and rotation in the auto-correlation function. It is shown that this 
representation is valid when N 2 «TB where N is the number of chips, T is the time duration of the 
signal, and B the bandwidth. 

In 1984 Costas introduced what came to be known as Costas Codes, codes whose hit 
array has a maximum side lobe level of 1 [57]. A good overview of different methods of 
constructing Costas Codes has been written by Moreno, Games and Taylor [58]. Another such 
paper was written by Golomb and Taylor [50, 59]. Variations on Costas Codes have been 
suggested, such as “Costas Signals with Intrapulse Frequency Modulation” [60]. In this paper, 
each individual chip of the code is frequency modulated. The peak to side-lobe ratio is the same, 
but the main lobe is narrower. Ricker has shown that if hops are allocated logarithmically in 
frequency, instead of linearly, then the autoambiguity function is the same regardless of center 
frequency and bandwidth [61]. If hops are allocated linearly then the sidelobes interfere with each 
other in different ways depending on the center frequency and bandwidth. 

One criteria often used to optimize signals is their cross-correlation properties. Freedman 
and Levanon have shown that “Any Two NxN Costas Signals Must Have at Least One Common 
Ambiguity Sidelobe if N > 3" [62]. Chang and Scarbrough have described a method of obtaining 
Costas arrays with a minimum number of cross-coincidences for a limited Doppler shift [63]. For 
an N x N Costas code, then q Costas codes can be found with a maximum cross-coincidence of 
one for a maximum Doppler shift of M where: (M+l)q<=N. Titlebaum and Marie have shown 
that if two cross coincidences are allowed, then Costas codes with minimum(2) cross 
coincidences can be found without limiting the allowed Doppler shift [64]. Certain constructions 
of Costas codes are shown to be better than others in this regard. In a separate paper they have 
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shown that certain sized Costas codes have better cross-correlation properties [65]. More work 
on this topic has been done by Drumheller and Titlebaum [66]. 

Another approach to optimizing the shape of the autoambiguity function is what is known 
as the Twin Processor or Uncertainty Product Function. In this approach more than one signal is 
transmitted into the environment, the scattering function is computed for each one, and then the 
two scattering functions are multiplied together. The effect is that real scatterers will add 
coherently, and spurious scatterers will be canceled. The drawback is that if the two signals are 
separated in time, then the velocity resolution is severely limited, and if they are separated in 
frequency the spacial resolution is severely limited. Work of this kind has been done by Mehta 
and Titlebaum [67], and Gustafson [69]. Naparst has worked on a system for combining 
information from several pings, where he uses the wideband processing and optimizes the shape 
of the autoambiguity function [69]. 

Signals have also been optimized to ignore, or “tolerate”, certain conditions. Altes has 
derived signals which are acceleration tolerant [70] and velocity tolerant [71]. Ricker has derived 
signals which are Doppler tolerant [72]. Altes has also designed signals for use when Doppler 
shifts are bounded [73]. 

Research has been done on echolocation signals used by bats, in hopes that evolution has 
in some way optimized the signals used in ways that are useful to sonar designers. Fenton and 
Bell found that an LFM followed by a constant frequency portion is correlated to an increased 
detection range in the species Myotis Volnas [74]. Kalko and Schmtzer found that certain bats 
switch from narrowband signals in free flight to wideband signals when approaching prey or 
cluttered areas (such as cave entrances) [75]. Saillant et al. have created a computational model 
of the cochlea and echolocation system of the bat [76]. They claim that this model has higher 
resolution than simple cross correlation. 
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APPENDIX B 

Alternatives to Computing the Probability of False Alarm 

In Chapter 3 it was shown that there are numerical problems computing the Heaviside 


coefficients of: 




Recall that 


A ■ 771 md A = * 


If all of the A k ‘s were precisely equal the corresponding probability function for the false alarm 
would be a Chi Squared distribution. It is now shown that the real characteristic function can be 
written as the Chi Squared distribution with a small correction. 


N J 

K ” n 1 -j2a>A kn 


where 


Ak » X k + l 1+1/A* 


if Xk is »1 then 1/7* «1. To emphasize this we will rename 1 / X k - £ k 


A k — 

1 + 6 y 
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The characteristic function now becomes 
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The beginning term (1 - jlcoy N is the characteristic function of a Chi Squared distribution of 
order N. The term it is multiplied by is one plus a small number, divided by one plus a small 
number, which is a number very close to one. This is highly suggestive of the fact that the false 
alarm may be approximated by a Chi Squared distribution. We will now show that 

|l/(l-y'2o>)|< 1. 

1 _ 1 + j2to 

1 - j2co 1 + 4 co 1 

(l + 4a> 2 ) l/2 
1+ 4<y 2 


Expanding the expression for <j>. 
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Taylor expanding yields 
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Dropping terms of order £ k 
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We now name these two terms: 


K = 0 Z 2 + A 


^ = (1 - j2a>y N 


H<t> lo = (l-j2a>y N 


j2a 
1 - j2co 



One approach we could take is to simply use the Chi Squared distribution, and drop any 
correction term. The ROC will still have information about the individual £ k ' s because the 
probability of detection can be computed exactly. 
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The probability of false alarm for a chi squared probability density has already been 
derived by Van Trees 30. 



l-/r 


r 


VM+1 


,M 


where I r is the incomplete gamma function which Van Trees defines as: 

uV M +1 M 

1 ^e-’dx 


It is important to note that this not the way that the incomplete gamma function is defined in 
Matlab. Van Trees also derives a approximation for large y . 

Fz 2 approx V Ml( 1- M/y"’) 


The notation with regards to y is a bit tricky. The likelihood ratio used in this paper (equation 


17) is equivalent to equation 398 on page 108 of Van Trees 30. Thus our y is equivalent to his 
y '. Van Trees then goes on to define: 


■„_ K + Q . 
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where <J 2 is the standard deviation of the noise and CT 2 is the standard deviation of the signal. In 
this work we have whitened the noise, so we take o\ = 1. The signal variance then corresponds 
to the eigenvalues of the whitened received signal <j* = A s , where the bar denotes the average. 
Preliminary numerical results have shown that the chi squared approximation is not only a good 
approximation, but an upper bound and thus the computed ROC will never better than the true 
ideal. No analytic method has yet been found to confirm this. A typical Probability of false alarm 
for six eigenvalues between 5 and 50 is shown in figure ? in comparison to the exact Chi Squared 
distribution, and the approximate Chi Squared distribution: 




Figure 31 Different expressions for the Probability of False Alarm 

We will now calculate the correction to the Chi Squared Probability of false alarm due to 
A <f> k . Defining s T = \[e k 
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The corresponding probability distribution is 
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making the substitution s = jco , ds - jdco , d(0 - ds / j 
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writing in terms of an inverse Laplace transform o£ 1 { } 
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This transform can be found in standard tables. It is: 
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The corresponding correction to the false alarm is 
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The correction to the false alarm is: 
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Alternatives to Computing the Probability of False Alarm: Advanced Methods of 
Computing the Heaviside Expansion 

A number of different methods of computing the Heaviside coefficients have been 
developed, unfortunately, the vast majority of them focus on computational efficiency and not on 
numerical accuracy. Also all of them compute coefficients which are equal to the coefficients 
found in the traditional manner, and are thus still limited by the machines precision. These 
schemes usually find some way to express the expansion in terms of a matrix inversion. 
Unfortunately, the case where all the roots are nearly equal correspond to a matrix whose rows 
are all nearly equal, which is also numerically difficult to compute. Mahoney and Sivazlian 76 
wrote a review of these methodologies as of 1983. Chang has developed a non-matrix recursive 
method 77, however the coefficients are still equivalent to the traditional coefficients. Matrix 
algorithms have been developed by Chang and Mott 78, and Uraz and N-Nagy 79. Deakin has 
written a paper on the history of the heaviside expansion 80. 
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